#include "mex.h"
#include <math.h>


static void MCPrice(double p[], double n1[], double n0[], int row, int c, double a, double b, double m, double sigma0, double x0, double t, double T, double r, double K)
{
        int i,j;
        double sigma1, dt=(T-t)/5000;
        double oucurrent, ouprevious;
        double xcurrent, xprevious;

        for (j=0; j<c; j++)
		{
			ouprevious = sigma0;
            xprevious = x0;

			for (i=0; i<row-1; i++)
            {
                oucurrent = ouprevious + a*(m-ouprevious)*dt+b*sqrt(dt)*n1[j*(row-1)+i];
                
                /*sigma1 = (2+sin(xprevious))*ouprevious;*/
                
                sigma1 = 0.3;
                
                xcurrent = xprevious + r*xprevious*dt+sigma1*xprevious*sqrt(dt)*n0[j*(row-1)+i];
                
                ouprevious = oucurrent; xprevious = xcurrent;
            }
            if(xcurrent >= K)
                p[j] = xcurrent - K;
            else
                p[j]=0;
		}
        


}



void mexFunction(int nlhs,mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
        double *p, *n0, *n1;
        int row,column;
        double a, b, m, sigma0, x0;
        double t, T, r, K;
        row=mxGetM(prhs[0])+1; column=mxGetN(prhs[0]);
        
        a=mxGetScalar(prhs[2]);
        b=mxGetScalar(prhs[3]);
        m=mxGetScalar(prhs[4]);
        sigma0=mxGetScalar(prhs[5]);
        x0 = mxGetScalar(prhs[6]);
        
        t=mxGetScalar(prhs[7]);
        T=mxGetScalar(prhs[8]);
        r=mxGetScalar(prhs[9]);
        K=mxGetScalar(prhs[10]);
        
        plhs[0]=mxCreateDoubleMatrix(1,column,mxREAL);
        
        p=mxGetPr(plhs[0]);
        n1=mxGetPr(prhs[0]);
        n0=mxGetPr(prhs[1]);

        MCPrice(p,n1,n0,row,column,a,b,m,sigma0,x0,t,T,r,K);


}
